Introduction
This project aimed to model brain blood vessel measurements of mice during a treatment protocol designed to elicit stress responses.
Measurement types included:
- Whisker stimulation response (vessel diameter before and after stimulation)
- Vessel centre and diameter pulsatility
- Red blood cell flow (i.e. speed and flux)
- Hypertensive challenge response (Correlation between blood pressure and diameter under different pressure conditions)
- Capillary density
- Vessel tortuosity
We believed that the mechanisms underlying each kind of measurement were independent, and moreover each measurement required different data filtering choices. We therefore conducted a separate analysis for each measurement type.
Our overall modelling approach broadly followed the recommendations in Gelman et al. (2020). For each analysis, we first constructed a simple mathematical model of the data generating process, then iteratively improved it, taking into account estimated predictive performance in and out of sample, simplicity, interpretability and computational feasibility.
We decided to employ a Bayesian modelling approach primarily because of the availability of non-experimental information, particularly structural information concerning the data making it important to partially pool information between categories given the relatively small number of measurements. In a Bayesian context partial pooling can be achieved using a prior distribution on the random effects. The capillary density dataset also benefitted from a Bayesian approach, as this made it straightforward to implement smoothing of adjacent capillary orders.
After verification of computational results, our analyses aimed to evaluate whether each of our models adequately captured the information contained in the available data, and then to interpret and present the results of the best models.
Overall strategy
Although our project involved several statistical analyses, we used a similar overall strategy in each case. This section describes the aspects of this strategy that were common to all of our analyses.
Data
This analysis involved three raw datasets, which we call the “ablation” dataset, the “hypertension” dataset and the “angio-architecture” dataset.
The ablation and hypertension datasets consist of real-valued measurements each with multiple categorical features, namely
- age of the measured mouse (adult or young)
- identity of the measured mouse
- stage of the treatment protocol when measured
- measured vessel type (penetrating artery, sphincter, bulb, first order capilary, etc)
The ablation and hypertension datasets contain different types of measurement. The ablation dataset includes measurements of whisker stimulation response, pulsatility and blood flow at various points of the treatment protocol. The hypertension dataset includes measurements of pressure to diameter correlation before and after two rounds of induced hypertension.
The angio-architecture dataset consists of real-valued measurements of vessel density and tortuosity.
Data processing
We ignored data from one mouse (id 310321) that was determined to be an outlier. 310321 is a mouse where we did not see any whisker response, it reacted to angiotensin II, but the BP increase was abrupted for a short while and then re-established. Perhaps due to a clot or a bubble in the venous catheter. This resulted in a biphasic and slow BP increase.
In all of our analyses we assumed that any missing measurements were caused by factors that were unrelated to our main target process, or in other words that the absent measurements were “missing at random”. We therefore did not attempt to model the measurement removal process explicitly.
Modelling approach
All of our models had a common structure, with generalised linear models used to describe information from measurements and multi-level prior distributions used to describe non-experimental information. The modelling choices open to us concerned the following questions:
What generalised linear model to use to model measurements?
Which interaction effects to estimate?
What structure to use for the prior model, and in particular how and to what extent to pool information between categories?
What quantitative values to use for our prior model?
In cases where a phenomenon of interest was assessed using multiple, potentially related measurements, whether to model the possible relatedness?
In order to answer these questions for each analysis, we started with a simple but plausible model, then iteratively added and removed components as described in Gelman et al. (2020). Our general aim was to achieve a better quantitative and qualitative description of the data generating process while avoiding computational issues. In particular, we focused on the estimated out of sample predictive performance as measured by the estimated leave-one-observation-out log predictive density (Vehtari, Gelman, and Gabry 2017) and qualitative agreement between predictive and observed observations in graphical checks.
In order to implement graphical predictive checking, we calculated quantiles of our models’ prior and posterior predictive distributions and plotted these alongside the measurements. We then inspected the graphs to ascertain whether the measurements were generally consistent with the predictions.
We calculated estimated leave-one-observation-out log predictive desities using the Python package arviz (Kumar et al. 2019).
Software
The raw data are found in csv files which are available from our project’s github repository at the following urls:
- https://github.com/teddygroves/sphincter/blob/main/data/raw/hyper_challenge.csv
- https://github.com/teddygroves/sphincter/blob/main/data/raw/data_sphincter_paper.csv
For each analysis, we conducted filtering and reshaping operations using the standard scientific Python stack and validated the resulting prepared datasets against custom data models constructed using the Python libraries pydantic (Pydantic developers 2022) and pandera (Niels Bantilan 2020). These models can be inspected at this url: https://github.com/teddygroves/sphincter/blob/main/sphincter/data_preparation.py.
Statistical computation was carried out using the probabilistic programming framework Stan (Carpenter et al. 2017) via the interface cmdstanpy (Stan Development Team 2022).
Analysis and serialisation of posterior samples was carried out using the Bayesian inference library arviz (Kumar et al. 2019).
Our analysis was orchestrated by the Python package bibat (Groves 2023).
Validation of statistical computation
We validated our statistical computation using standard Hamiltonian Monte Carlo diagnostics, including the improved \(\hat{R}\) statistic (Vehtari et al. 2021) as well as inspection for post-warmup divergent transitions (Betancourt 2017), problematic EBFI statistics, tree depth or effective sample size to total sample size ratios. All reported models had improved \(\hat{R}\) close to 1 and no divergent transitions or other signs of algorithm failure.
Reports of findings
In general the findings we report were cases where our best performing model’s posterior probability distribution gave a high probability to an interesting and interpretable quantity having an interesting value. Since most of the interpretable parameters in our models are normalised so that zero indicates no effect, in practice our main findings are mostly that a parameter, or the difference between some parameters, is likely to be substantially different from zero according to our best model’s posterior distribution. We present results with this form using histograms of the relevant marginal posterior distributions.
Reproducibility
See the repository readme for instructions on reproducing our analysis: https://github.com/teddygroves/sphincter/blob/main/README.md
Main findings
Whisker stimulation
Our analysis of whisker stimulation data indicates that the hypertension and sphincter ablation treatments are associated with lower whisker stimulation response, as measured by log diameter change, compared with the baseline treatment. The effect from whisker stimulation is greater than from hypertension.
Figure Figure 1 illustrates this finding by showing the distribution of posterior samples for treatment effects relative to baseline from our best whisker stimulation model.
Our analysis did not indicate any substantial difference between old and adult mice, or any noticeable vessel type:treatment interaction effects. This can be seen from figure Figure 2, which shows posterior quantiles for age and vessel type:treatment interaction effects in our model that included both of these.
We found some difference between vessel type effects: sphincters had the greatest relative diameter change in response to whisker stimulation, and bulbs the smallest. Figure Figure 3 shows these.
Pulsatility
Our analysis of vessel centre and diameter pulsatility yielded the following conclusions:
Adult mice have higher vessel diameter pulsatility than old mice, whereas old mice have slightly higher centre pulsatility.
Sphincter ablation correlates with increased diameter pulsatility, with no strong interaction effects. On the other hand there is no clear effect of sphincter ablation on centre pulsatility.
Figure 4 plots the distribution of age effect differences (adult minus old) for each measurement type in our final model. This graph shows that, in this model, the age effect for adult mice was higher than for old mice in every single posterior sample: in other words, according to our model older mice have lower diameter pulsatility. There is a smaller opposite trend for centre pulsatility measurements, but it is not clearly separated from zero, indicating that the direction of the effect is not fully settled.
Figure 22 shows the distribution of posterior draws for sphincter ablation effects compared with the immediately prior protocol stage (“after hypertension”). The ablation/diameter parameter is greater than the after hypertension/diameter parameter in 98% of posterior samples, whereas there is no clear effect on centre pulsatility.
Red blood cell flow
Our main result regarding red blood cell flow is that both RBC speed and flux were higher in the adult mice compared with the old mice. Figure Figure 6 illustrates this finding by plotting the relevant marginal posterior histograms.
We also quantified treatment effects on red blood cell flow, as shown in Figure 7:
Hypertensive challenge
Our hypertensive challenge data also showed pronounced age and treatment differences, as shown in figure Figure 8. Specifically, we found that, for the adult mice, blood pressure and vessel diameter were less correlated, and that the hyper2 treatment was associated with increased pressure-diameter correlation compared with the hyper1 treatment.
Capillary density
The capillary density dataset had a somewhat different structure to the other data, with no treatments, more vessel types, with clear correlation between measurements corresponding to adjacent vessel types. We therefore used a different statistical approach, with smoothing components for parameters of adjacent vessel types. See Section 8.2 for more about this model.
Our analysis indicated that the old mice tended to have lower density than the adult mice for capillaries of order 1 to 5 and higher density for capillaries of order 10 to 12, and for ascending venules, as shown in figure Figure 9.
Vessel tortuosity
We also fit a smoothed regression model to measurements of vessel tortuosity from the angio-architecture dataset. See Section 9.2 for details about this model.
As shown in Figure 10, our analysis showed a clear difference between the adult and old mice in the upper vasculature, with the old mice tending to have substantially more tortuous pial arteries and penetrating arterioles. Other vessels were similarly tortuous for both age categories.
Details of the whisker stimulation analysis
In order to measure how vascular responsiveness changed during our experimental protocol, the diameters of different vessel types were recorded before and during whisker stimulation, at baseline, post-hypertension and post-ablation stages. We aimed to assess statistical relationships between the known factors and the relative change in vessel diameter before and after stimulation.
In particular, we were interested in differences between old and adult mice and in the effect of sphincter ablation.
Dependent variable
The ratio of the peak compared with the pre-stimulation level for each mouse at each stage, on natural logarithmic scale, also known as the ‘log change’, was standardised by subtracting the overall mean and dividing by the standard deviation, then treated as a single measurement. This way of the measurements was chosen to facilitate modelling, as log change is a symmetric and additive measure of relative change (see Tornqvist, Vartia, and Vartia (1985)).
Note that when the difference between the two values \(v1\) and \(v2\) is far less than 1, the log change \(\ln{\frac{v2}{v1}}\) is approximately the same as the more widely used relative difference measure \(\frac{v2-v1}{v1}\).
Statistical Models
The final model is shown below:
\[\begin{align} y_{vtm} &\sim ST(\nu, \hat{y}_{vtm}, \sigma_{t}) \label{eq-whisker-model} \\ \hat{y}_{vtm} &= \mu_a \nonumber \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_v \nonumber \\ \alpha^{treatment}_t &\sim N(0, \tau^{treatment}) \nonumber \\ \alpha^{vesseltype}_v &\sim N(0, \tau^{vesseltype}) \nonumber \\ \nu &\sim Gamma(2, 0.1) \nonumber \\ \sigma_t &\sim HN(0, 0.5) \nonumber \\ \mu &\sim N(0, 0.7) \nonumber \\ \tau^{treatment} &\sim HN(0, 0.7) \nonumber \\ \tau^{vesseltype} &\sim HN(0, 0.7) \nonumber \end{align}\]
In equation \(\eqref{eq-whisker-model}\), the term \(ST\) indicates the student t distribution, \(N\) indicates the normal distribution, \(Gamma\) the gamma distribution and \(HN\) the ‘half-normal’ distribution, i.e. the normal distribution with support only for non-negative numbers. Subscripts indicate indexes and superscripts indicate parameter labels.
This model has independent effects for treatments (\(\alpha^{treatment}\)), vessel types (\(\alpha^{vessel}\)) and age (\(\mu\)). The treatment and vessel type effects have hierarchical priors, reflecting the need to partially pool information between different treatment and vessel type categories. This structure allows our model to accommodate the full spectrum between different categories being highly heterogenous—in this case the corresponding \(\tau\) parameter will be large—and on the other hand high similarity, leading to small \(\tau\) parameters. The student t distribution was chosen to reflect the heavy tails we observed in the data, with the parameter \(\nu\) controlling the level of heaviness.
The prior standard deviation 0.7 was chosen because it led to what we judged to be a reasonable allocation of prior probability mass over possible data realisations. The prior for the student t degrees of freedom parameter \(\nu\) was set following the recommendation in Juárez and Steel (2010).
As well as this model, we also present results from a more complex model that adds vessel type:treatment and age:treatment interaction effects. The full specification of this “big” model is as follows, using the same notation as equation \(\eqref{eq-whisker-model}\):
\[\begin{align} y_{vtm} &\sim ST(\nu, \hat{y}_{vtm}, \sigma_{t}) \label{eq-whisker-model-interaction} \\ \hat{y}_{vtm} &= \mu_a \nonumber \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_v \nonumber \\ &+ \alpha^{vesseltype:treatment}_{vt} \nonumber \\ &+ \alpha^{age:treatment}_{vt} \nonumber \\ \alpha^{treatment}_t &\sim N(0, \tau^{treatment}) \nonumber \\ \alpha^{vessel}_v &\sim N(0, \tau^{vessel}) \nonumber \\ \alpha^{vesseltype:treatment}_{vt} &\sim N(0, \tau^{vesseltype:treatment}) \nonumber \\ \alpha^{age:treatment}_{at} &\sim N(0, \tau^{age:treatment}) \nonumber \\ \nu &\sim Gamma(2, 0.1) \nonumber \\ \sigma &\sim HN(0, 0.5) \nonumber \\ \mu &\sim N(0, 0.7) \nonumber \\ \tau^{treatment} &\sim HN(0, 0.7) \nonumber \\ \tau^{vessel} &\sim HN(0, 0.7) \nonumber \\ \tau^{age:treatment} &\sim HN(0, 0.7) \nonumber \\ \tau^{vesseltype:treatment} &\sim HN(0, 0.7) \nonumber \end{align}\]
Results
Figure 11 shows the observed log change measurements with colours illustrating the various categorical information. Note that the measurements have different dispersion depending on the treatment, indicating a the need for a model with distributional effects.
Figure 12 compares the measurements with our model’s posterior predictive distribution. This shows that our model achieved a reasonable fit to the observed data. There is a pattern in the model’s bad predictions, in that these tend to be for higher baseline measurements. However, we judged that this pattern was small enough that for our purposes we could disregard it.
To back up our finding that there were no important interaction effects, we compared the predictive performance of our final model whisker-ind with whisker-big, the best-performing model with more interactions. The results are shown below in figure Figure 13:
whisker-ind and the best performing interaction model whisker-big.
The two models have similar estimated predictive performance, indicating that there is no gain from considering interaction effects in this case.
Details of the pulsatility analysis
The pulsatility data consisted of fast measurements of diameter and center point for the same mice. These measurements were Fourier-transformed, and the harmonics of the transformed data were interpreted as representing the pulsatility of the measured quantities.
Dependent variable
We used the first harmonic of each transformed time series as a dependent variable. It might have been preferable to aggregate all the available power harmonics, but this would have complicated our measurement model, and in any case power at the subsequent harmonics was typically negligible compared with the first.
Questions
As well as the results reported in the main findings section, we were also interested in these additional questions:
Question a How does blood pressure affect diameter and centre pulsatility?
Question b Do hypertension and sphincter ablation influence diameter and centre pulsatility differently?
Description of the dataset
As well as the categorical data described above, our pulsatility analysis also took into account measurements of each mouse’s blood pressure at the femoral artery.
The final dataset included 514 joint measurements of diameter and centre pulsatility, calculated as described above. These measurements are shown in Figure 14.
Figure 15 shows the relationship between pressure and the measurements in our dataset for both age categories. The light dots show raw measurements and the darker dots show averages within evenly sized bins.
Figure 16 shows the relationship between diameter and the measurements in our dataset for all vessel type categories. The light dots show raw measurements and the darker dots show averages within evenly sized bins. There is a clear positive relationship between measured absolute diameter and diameter pulsatility, and it is approximately the same for all vessel types.
Statistical models
We knew from prior studies that the power harmonics should individually follow exponential distributions [REFERENCE FOR THIS]. This consideration motivated the use of exponential generalised linear models for both the centre and diameter pulsatility measurements. In this model, given measurement \(y\) and linear predictor \(\eta\) the measurement probability density is given by this equation:
\[\begin{align} p(y\mid\eta) &= Exponential(y, \lambda) \label{eq:pulsatility-measurement-model} \\ &= \lambda e^{-\lambda y} \nonumber \\ \ln{\frac{1}{\lambda}} &= \eta \label{eq:link-function} \end{align}\]
The log link function \(\eqref{eq:link-function}\) was chosen so that linear changes in the term \(\eta\) induce multiplicative changes in the mean \(\frac{1} {\lambda}\) of the measurement distribution, as we believed the effects we wanted to model would be multiplicative.
We compared four different ways of parameterising \(\eta\) based on the information available about a given measurement, corresponding to three hypotheses about the way the data were generated.
The simplest model, which we labelled “basic”, calculates the linear predictor \(\eta^{basic}_{vtad}\) for an observation with vessel type \(v\), treatment \(t\), age \(a\) and diameter \(d\) as follows:
\[\begin{align} \label{eq:basic} \eta^{basic}_{vtad} &= \mu_{a} \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{v} \nonumber \\ &+ \beta^{diameter} \cdot \ln{d} \nonumber \end{align}\]
The basic model provided a plausible baseline against which to compare the other models.
Next we constructed a more complex model by extending the basic model with interaction effects, resulting in the following linear predictor:
\[\begin{align} \label{eq:interaction} \eta^{interaction}_{vtad} &= \mu_{a} \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{d,vesseltype(n)} \\ &+ \alpha^{treatment:vesseltype}_{tv} \nonumber \\ &+ \beta^{diameter}_{d} \cdot \ln{d} \nonumber \end{align}\]
Next we constructed a model that adds to the basic model parameters that aim to capture possible effects corresponding to the blood pressure measurements. To compensate for collinearity between age and pressure, our “pressure” model does not use the observed pressure as a predictor, but rather the age-normalised pressure, calculated by subtracting the mean for each age category from the observed pressure measurement. The model for the linear predictors $ ^{pressure}_{vatdp}$ with age-normalised pressure measurement \(p\) is then
\[\begin{align} \label{eq:pressure-model} \eta^{pressure}_{vatdp} &= \mu_{a} \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{v} \nonumber \\ &+ \beta^{diameter}_{d} \cdot \ln{d} \nonumber \\ &+ \beta^{pressure}_{a} \cdot p \nonumber \end{align}\]
Finally, we made a model that includes a pressure effect but no age-specific parameters from the pressure model. This was to test whether any age effects are due to the collinearity between age and pressure. The pressure-no-age model’s linear predictors \(\eta^{pressure\ no\ age}_{vatdp}\) are calculated as shown in equation \(\eqref{eq:pressure-no-age-model}\). Note that, unlike in equation \(\eqref{eq:pressure-model}\), the \(\mu\) and \(\beta^{pressure}\) parameters in equation \(\eqref{eq:pressure-no-age-model}\) have no age indexes.
\[\begin{align} \label{eq:pressure-no-age-model} \eta^{pressure\ no\ age}_{vatdp} &= \mu \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{v} \nonumber \\ &+ \beta^{diameter}_{d} \cdot \ln{d} \nonumber \\ &+ \beta^{pressure} \cdot p \nonumber \end{align}\]
In all of our models the \(\alpha\) parameters were given independent, semi-informative, hierarchical prior distributions to allow for appropriate information sharing. The \(\beta\) and \(\mu\) parameters were given independent, semi-informative, non-hierarchical prior distributions.
Results
We estimated the leave-one-out log predictive density for each model using the method described in Vehtari, Gelman, and Gabry (2017) and implemented in Kumar et al. (2019). The results of the comparison are shown below in Figure 17.
We evaluated our models’ fit to data using prior and posterior predictive checking, with the results for the pressure model shown in Figure 18.
Inspecting of the interaction model output showed that none of the interaction effect parameters that differed substantially from zero, as can be seen in Figure 19.
From this result, together with the worse estimated out of sample predictive performance as shown in Figure 17, we concluded that there were no important interaction effects, so that we could essentially discard the interaction model.
Figure 20 shows the marginal posterior distributions for other effect parameters in all three models. Note that the parameters b_diameter are strongly positive for diameter pulsatility in all models and also mostly positive for centre pulsatility. There is also a strong trend for diameter pulsatility to decrease with the order of the vessel and no particular vessel type trend for centre pulsatility.
Answers to specific questions
The poorer estimated out of sample predictive performance of the pressure-no-age model compared with the other models, as shown in Figure 17, indicates that our pressure measurements did not fully explain the observed difference between adult and old mice. It is nonetheless possible that different pressure explains the difference between old and adult mice, but that the pressure measurements did not reflect the true pressure at the measured vessels. This is plausible since the pressure measurements were taken at a different location.
Figure 21 shows the difference in \(\beta^{pressure}\) parameters for old and adult mice in the pressure model in order to answer Question a. This shows a weak tendency of the pressure effect on diameter pulsatility to be more positive for adult mice than for old mice, and a strong opposite tendency for centre pulsatility. Taking the absolute values into account, the analysis suggests that greater measured pressure is not strongly related to diameter pulsatility and correlates with reduced centre pulsatility for adult mice but not for old mice.
To illustrate the effect of treatments, and specifically sphincter ablation relative to hypertension (i.e. to answer Question b) Figure 22 shows the difference between the effect for each treatment and the baseline treatment effect. There is a clear effect of ablation to increase diameter pulsatility and no clear effects of hypertension on diameter pulsatility or of either treatment on centre pulsatility.
To get an idea about how the effect of sphincter ablation on diameter pulsatility compares quantitatively with the effect of hypertension, we fit the basic model to the full dataset, without excluding measurements from either hypertension treatment. Figure 23 shows the main result from fitting this model: ablation and hypertension had similarly positive effects on diameter pulsatility. Interestingly there is no clear effect from the second hypertension treatment.
Details of the red blood cell flow analysis
Our measurements included flow data recording the measured speed and flux of red blood cells through some vessels. This data is interesting because it allows for inference of the local blood pressure, which determines the speed.
We were interested in whether the speed and flux tended to be different between old and adult mice for a given vessel type and treatment, as this would indicate that the pressure would likely be similar as well.
Data processing
There is a significant missing data issue in this case: both speed and flux measurements were available for only 271 out of 1525 raw measurements. We therefore conducted separate analyses for speed and flux even though we suspected that these two quantities are related.
Dependent variable
We modelled speed and flux measurements on natural logarithmic scale as we expected multiplicative effects and this transformation ensures support on the whole real number line, simplifying modelling.
The resulting measurements are shown in figure Figure 24.
From a glance at these graphs it is clear that there were treatment effects on both measurement types, particularly from the hypertension treatments, that there are treatment-related distributional effects, and that both speed and flux reduce as vessel order increases.
Statistical models
As in the whisker case we investigated the results of fitting models with and without interaction effects. Again we found no large or fully resolved interactions and therefore used the smaller model for further analysis.
Our final model had this specification:
\[\begin{align} \ln(y_{vtm}) &\sim N(\hat{y}_{vtm}, \sigma_t) \label{eq-flow-model} \\ \hat{y}_{vtm} &= \mu_a \nonumber \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{v} \nonumber \\ \alpha^{treatment}_t &\sim N(0, \tau^{treatment}) \nonumber \\ \alpha^{vesseltype}_v &\sim N(0, \tau^{vesseltype}) \nonumber \\ \sigma_t &\sim HN(0, 0.5) \nonumber \\ \mu &\sim N(0, 0.5) \nonumber \\ \tau^{treatment} &\sim HN(0, 0.5) \nonumber \\ \tau^{vesseltype} &\sim HN(0, 0.5) \nonumber \end{align}\]
We chose the prior standard deviation 0.5 after prior predictive checking to ensure reasonably tight coverage of the observed data.
For investigation of interaction effects we fit another model that extended our final model with a vessel type:treatment interaction effect as follows:
\[\begin{align} \hat{y}_{vtm} &= \mu_a \nonumber \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{v} \nonumber \\ &+ \alpha^{vesseltype:treatment}_{vt} \nonumber \\ \alpha^{vesseltype:treatment}_v &\sim N(0, \tau^{vesseltype:treatment}) \nonumber \\ \tau^{vesseltype:treatment} &\sim HN(0, 0.5) \nonumber \end{align}\]
Results
Our statistical model successfully captured all of these trends, as can be seen from figure Figure 25
Interactions
To test for important interaction effects we fit a model with vessel type:treatment interaction parameters. This model achieved marginally worse loo elpd scores as shown in figure Figure 27.
For completeness the interaction effects are shown in figure Figure 28. No effects are clearly separated from zero. The sphincter/ablation effect on RBC speed is notably different from the others. We fit several models with sparsity-inducing priors including the regularised horseshoe Piironen and Vehtari (2017) to see if it was possible to resolve this effect, but were unsuccessful. From this we conclude that any real effect is too small to be easily detected in our dataset.
Details of the hypertensive challenge analysis
Dependent variable
The raw data on hypertension challenge were correlation coefficients relating blood pressure and vessel diameter, which are constrained to lie on the \([-1, 1]\) interval. For easier modelling we transformed these by applying an inverse hyperbolic tangent function for use in modelling. The dependent variables then had support on the entire real number line.
The resulting dataset is shown in figure Figure 29. The transformed correlation coefficients do not appear extremely dispersed, indicating that standard modelling techniques ought to be able to describe them.
Statistical model
Our final statistical model had the following form.
\[\begin{align} \ln(y_{vtm}) &\sim N(\hat{y}_{vtm}, \sigma_v) \label{eq-hypertension-model} \\ \hat{y}_{vtm} &= \mu_a \nonumber \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{v} \nonumber \\ \alpha^{treatment}_t &\sim N(0, 0.5) \nonumber \\ \alpha^{vesseltype}_v &\sim N(0, \tau^{vesseltype}) \nonumber \\ \sigma_v &\sim HN(0, 0.5) \nonumber \\ \mu &\sim N(0, 0.5) \nonumber \\ \tau^{vesseltype} &\sim HN(0, 0.5) \nonumber \end{align}\]
This model is different from the others in that we did not partially pool the treatment effects, since there were only two of these in this case. We also allowed the measurement error parameters \(\sigma\) to vary according to vessel type, since this improved model fit and predictive performance.
As in the other analyses, for investigation of interaction effects we fit another model that extended our final model with a vessel type:treatment interaction effect as follows:
\[\begin{align} \hat{y}_{vtm} &= \mu_a \nonumber \\ &+ \alpha^{treatment}_{t} \nonumber \\ &+ \alpha^{vesseltype}_{v} \nonumber \\ &+ \alpha^{vesseltype:treatment}_{vt} \nonumber \\ \alpha^{vesseltype:treatment}_v &\sim N(0, \tau^{vesseltype:treatment}) \nonumber \\ \tau^{vesseltype:treatment} &\sim HN(0, 0.5) \nonumber \end{align}\]
Results
Figure 30 shows that, as in the other cases, including interaction effects did not improve estimated predictive performance. The expected leave-one-observation-out predictive density is estimated to be higher for the model hypertension-basic, which did not have any interaction effects, than for the interaction model hypertension-big. Furthermore, the standard error for the difference in these two quantities is much smaller than the difference itself, indicating that the ranking is robust.
hypertension-big model is clearly worse.
Figure 31 shows the marginal distributions for the non-hierarchical parameters in our final model. As shown in Figure 8, there are clear differences between the different mu and a_treatment parameters, indicating that there are important age and treatment effects. There is also a difference between parameters corresponding to penetrating arterioles and the other vessel type parameters. This reflects the fact that the penetrating arterioles had lower pressure vs diameter correlation coefficients than other vessels, and that the penetrating arteriole measurements were more dispersed.
Figure 32 shows graphical prior and posterior predictive checks for our final hypertension model. The fit is fairly good, with no obvious systematic pattern in the errors, though slightly more observations lie outside the plotted intervals than might be expected.
Details of the density analysis
The density dataset consisted of 144 measurements from five adult and four old mice.
Dependent variable
The dependent variable in this case was capillary density, measured as in length of vessel per unit of volume. These measurements are shown in Figure 33 on both natural and logarithmic scales.
We noticed several interesting patterns in this data.
First, there is a clear trend for capillary density to increase with vessel order until the cap6 vessel type and then to decrease, with adjacent vessel types tending to have similar densities.
Second, we noticed that the measurements are somewhat more dispersed for higher-order capillaries, with the level of dispersion increasing somewhat smoothly.
Statistical model
We modelled the measurement process using a linear model on logarithmic scale. In this model, given measurement \(y\), linear predictor \(\hat{y}\) and a vessel type specific measurement error parameter \(\sigma\), the measurement probability density is given by this equation:
\[\begin{equation} p(y\mid\alpha, \sigma) = Normal(\ln{y}\mid \alpha, \sigma) \label{eq:measurement-model-density} \end{equation}\]
In order to capture the observed smoothness between adjacent vessel types, we used Gaussian random walk priors for the \(\alpha^{age:vesseltype}\) parameters:
\[\begin{align} \alpha_{a,v} &\sim N(\alpha_{a,v-1}, \lambda^{\alpha}_a)\label{eq:smooth} \\ \end{align}\]
We also used a Gaussian random walk prior on log standardised scale for the measurement error parameter \(\sigma\), where \(\sigma^{std}=\frac{\sigma}{sd(\ln y)}\):
\[\begin{align} \ln\sigma^{std}_{v} &\sim N(\ln\sigma^{std}_{v-1}, \lambda^s)\label{eq:smooth-sigma} \\ \end{align}\]
This approach to smoothing parameters corresponding to ordered categories is essentially the same as that used in Gao et al. (2019) to model age effects on voting behaviour. As explained in that paper, the random walk priors allow for information sharing between categories, without the need for detailed assumptions about the functional form of the overall relationship.
The other priors in our model were as follows (units are on standardised logarithmic scale):
\[\begin{align} \alpha_{1} &\sim N(0, 0.1) \label{eq:density-other-priors} \\ \ln\sigma^{std}_{1} &\sim N(-1.5, 1) \nonumber \\ \lambda^{\alpha} &\sim N(0, 0.3) \nonumber \\ \lambda^{\sigma} &\sim N(0, 0.3) \nonumber \\ \end{align}\]
Results
Figure 34 shows the overall fit of our model to the observed data. We judged that the overall fit to the observed measurements was adequate and did not attempt model evaluation using estimated leave-one-out density as given the highly correlated data it would be difficult to perform the necessary estimates with acceptable accuracy. In particular, our model captured the increase and decrease in density with increasing capillary order for both age categories, and the higher dispersion at higher capillary orders.
Figure 35 shows posterior 94% credible intervals for the quantity \(\alpha_{2} - \alpha_{1}\) for each vessel type; in other words, the overall difference, on standardised log scale, in densities between old and adult mice for each vessel type.
There is a clear trend, with four low order capillary types separated from zero on the left and two high order types separated on the right.
Details of the tortuosity analysis
Like the density analysis, the tortuosity analysis modelled the architecture dataset.
Dependent variable
The dependent variable in this case was vessel tortuosity, measured as total vessel length divided by straight-line length.
The raw measurements are shown in Figure 36.
There is a clear outlier for mouse 270220 at the 4th capillary: we excluded this measurement but kept other measurements from this mouse as they did not seem particularly anomalous.
Otherwise, we observed a similar pattern of smoothness between adjacent vessel as with the density measurements, motivating the use of a model with smoothed latent parameters. Unlike with the density measurements, however, there is some variety in the size of the jumps from vessel type to vessel type. In particular, the adult mice show a sharp upward jump from vessel pa to the first capillary, and both the adult and old mice show a downward jump in tortuosity at the ascending venule.
Statistical model
As with the density dataset, we used a regression model with random-walk priors to smooth adjacent vessel type effects within age categories and measurement errors for adjacent vessel types. In order to appropriately capture the heterogeneity in jump sizes, we used a student-t distribution rather than a normal distribution for the random walk in vessel type effects, with the degrees of freedom parameter modelled hierarchically.
The model specification is as follows, with \(\ln{y}^{std}\) representing the log-transformed and then standardised tortuosity measurements:
\[\begin{align} p(y\mid\alpha, \sigma) &= Normal(\ln{y}^{std}\mid \alpha, \sigma) \label{eq:measurement-model-tortuosity} \\ \alpha_{a,v} &\sim ST(\nu, \alpha_{a,v-1}, \lambda^{\alpha}_a) \nonumber \\ \ln\sigma^{std}_{v} &\sim N(\ln\sigma^{std}_{v-1}, \lambda^s) \nonumber \\ \nu &\sim \Gamma(2, 0.1) \nonumber \\ \alpha_{,0} &\sim N(0, 1) \nonumber \\ \ln\sigma^{std}_{0} &\sim N(-1.5, 0.8) \nonumber \\ \lambda &\sim N(0, 0.2) \nonumber \\ \lambda^{s} &\sim N(0, 0.2) \nonumber \end{align}\]
Results
Figure 37 shows this model’s posterior predictive distribution alongside the retained measurements, indicating a reasonably good fit.
The posterior predictive distributions for upper vessels “pial artery” and “pa” appear different for old and adult mice. To ascertain the extent and certainty of this contrast we plotted the differences between vessel type effects for old and adult mice in Figure 38. From this plot it is clear that, according to our model, there is a pronounced and robust difference, with over 97% probability of each effect being greater for old mice.